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EXECUTIVE  SUMMARY 


We  have  demonstrated  robust  algorithms  for  generating  fits  in  Optical  Measure  Theory 
(OMT)  to  radiance  data  from  the  TIROS  Operational  Vertical  Sounder  (TOVS).  These 
algorithms  generated  ph > sically  meaningful  fits  on  100%  of  a  test  set  of  45  TOVS  radiance 
scans,  and  were  successful  both  on  short  wavelength  (700  cm'1)  and  on  long  wavelength  (2250 
cm1)  TOVS  data.  The  resulting  OMT  temperature  profiles  exhibit  meteorological 
characteristics  and  appear  to  be  suitable  for  determination  of  atmospheric  structure  parameters 
characterizing  the  large-scale  vertical  temperature  structure  of  the  atmosphere,  and  also  suitable 
for  generation  of  input  data  for  numerical  weather  prediction  codes. 
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1.  INTRODUCTION 


The  mathematical  difficulties  encountered  in  attempting  to  invert  satellite  radiance  data  to 
obtain  atmospheric  temperature  profiles  have  been  of  considerable  theoretical  interest  to  the 
mathematical  and  satellite  meteorology  communities.  In  this  report  we  outline  the  application  of 
a  new  approach  to  the  solution  of  this  problem,  Optical  Measure  Theory  (OMT),  which  has  been 
invented  by  Dr.  J.  I.  F.  King  of  the  Air  Force  Geophysics  Laboratory  (King  1989  ;  Leon  and 
King  1988). 

The  upwelling  radiance  is  determined  by  the  equation  of  radiative  transfer,  which  is  a 
Fredholm  integral  equation  of  the  first  kind,  integrating  over  the  Planck  function  as  a  function 
of  height.  The  mathematical  difficulties  in  carrying  out  inversions  stem  from  this  simple 
physical  description.  Since  integration  is  intrinsically  a  smoothing  operation,  the  radiance 
profile  will  tend  to  be  smoother,  qualitatively  speaking,  than  the  atmospheric  temperature 
profile  which  generates  it.  Furthermore,  since  the  solutions  of  these  integral  equations  are  not 
unique,  many  atmospheric  temperature  profiles  (in  fact,  formally  an  uncountably  infinite  set) 
can  give  rise  to  the  same  radiance  profile,  and  thus  cannot  be  distinguished  on  the  basis  of 
radiance  observations. 

Clearly,  a  radiance  profile  has  only  a  limited  information  content,  relative  to  the  problem  of 
reconstructing  the  atmospheric  temperature  profile.  (We  intend  to  make  this  statement 
mathematically  rigorous  in  later  research  by  an  information  theoretical  analysis,  which  will 
address  the  actual  information  content,  in  bits,  of  radiance  observations.)  Therefore, 
construction  of  an  atmospheric  temperature  profile  which  replicates  the  detailed  structure  of  an 
in  situ  radiosonde  observation  is  probably  unrealistic,  however,  extracting  information  about  the 
overall  morphology  of  the  atmospheric  temperature  profile  is  a  meaningful  objective  within  the 
context  of  present  theoretical  developments.  We  expect  that  this  overall  morphological 

information  concerning  the  atmospheric  temperature  profile  may  be  expressed  in  terms  of 
“atmospheric  structure  parameters”,  such  as  temperature  lapse  rates  and  tropopause  temperatures 
and  pressures,  which  are  meaningful  quantities  in  meteorological  research.  Furthermore,  these 
quantities  are  appropriate  input  parameters  for  numerical  weather  prediction  codes,  which  are  in 
any  event  insensitive  to  the  small-scale  atmospheric  temperature  structures  which  are  found  in 
radiosonde  temperature  profiles. 

OMT  has  a  number  of  desirable  theoretical  features  that  make  it  especially  appropriate  as  an 
approach  to  the  inversion  problem.  In  OMT,  a  smooth  functional  form  is  fit  to  the  radiance 
profile;  thus,  OMT  represents  one  approach  for  solution  of  the  inversion  problem  by 

regularization,  in  the  sense  that  the  space  of  possible  solutions  is  restricted  by  the  imposition  of 
reasonable  physical  constraints.  This  restriction  of  the  class  of  functions  used  in  representing 
the  radiance  profiles  is  also  reflected  in  a  restriction  of  the  class  of  functions  used  in  the 
representation  of  the  atmospheric  temperature  profiles.  In  effect,  we  select  the  smoothest 
possible  representation  of  the  radiance  profile  consistent  with  the  radiance  data  set  at  an 
appropriate  level  of  confidence.  Crudely  speaking,  this  represents  a  “least  information” 
representation  in  the  sense  that  we  have  imposed  a  minimum  of  information  not  actually  present 
in  the  radiance  measurements  in  the  process  of  carrying  out  our  smoothing  through  a  functional 
fit  to  the  radiance  data.  The  resulting  atmospheric  temperature  is  then  also  "bias-free"  in  the 
sense  that  it  contains  no  a  priori  information  concerning  the  temperature  structure  of  the 
atmosphere.  Such  a  representation  may  fairly  be  said  to  have  been  constructed  if  the  functions 

from  which  the  radiance  representation  is  constructed,  which  we  term  “basis  functions”,  are 

chosen  on  grounds  of  fundamental  physical  principles  describing  the  problem.  In  fact,  in  OMT 
this  choice  is  made  on  the  basis  of  the  mathematical  description  of  a  radiative  atmosphere. 
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It  is  important  to  note  that  in  OMT  we  have  imposed  the  constraints  of  functional 
smoothness  on  the  radiance  data,  that  is,  as  close  to  the  instrument  as  possible.  All  remaining 
transformations  in  OMT  required  for  the  calculation  of  atmospheric  temperature  profiles  are 
strictly  one-to-one  transformations  (and  hence  information-conserving).  OMT  is  theoretically 
advantageous  in  this  regard  also  because  the  relation  between  the  theoretical  results  of  OMT  and 
instrumental  characteristics  is  particularly  clear  and  straightforward  in  principle. 

The  functional  representation  for  the  atmospheric  temperature  profiles  is  calculated  from  the 
radiance  profile  basis  functions  by  a  generalization  of  the  Laplace  transform.  This 
transformation  is  motivated  by  realizing  that  for  a  large  class  of  atmospheric  weight  functions, 
the  equation  of  radiative  transfer  may  be  viewed  as  an  integral  transform,  specifically  a 
generalization  of  the  Laplace  transform.  This  observation  is  one  of  the  crucial  theoretical 
insights  contained  in  OMT.  As  a  result,  the  choice  of  basis  function'  used  in  the  representation 
of  the  radiance  profile  and  in  the  representation  of  the  atmospheric  temperature  profile  may  not 
be  made  independently;  the  temperature  profile  basis  functions  are  (generalized)  inverse  Laplace 
transforms  of  the  radiance  profile  basis  functions.  In  this  sense,  OMT  may  be  said  to  be  an 
algebra  of  radiance  profiles,  since  profile-like  functional  forms  are  fundamental  entities  of 
OMT. 

We  shall  demonstrate  in  this  report  how  OMT  may  be  used  to  construct  atmospheric 
temperature  profiles  which  have  an  overall  morphology  characteristic  of  observed  atmospheric 
temperature  profiles  and  which  may  be  sued  to  obtain  atmospheric  structure  parameters  suitable 
for  input  to  numerical  weather  prediction  codes.  In  the  study  reported  on  here,  the  radiance 
data  analyzed  were  obtained  from  the  TIROS  Operational  Vertical  Sounder  (TOVS),  however, 
the  general  principles  may  be  applied  to  other  atmospheric  sounders,  such  as  on  the  DMSP 
satellites.  We  have  demonstrated  robust  algorithms  which  yield  valid  OMT  solutions  from  TOVS 
radiance  data  with  reasonably  high  probability.  Additional  research  will  improve  the 
computational  efficiency  and  robustness  of  OMT  and  generate  further  applications  in  satellite 
meteorology  and  numerical  weather  prediction. 
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2.  THEORY 

2.1  Introduction  and  Development  of  Optical  Measure  Theory.  In  this  subsection,  we  include 
for  completeness  a  derivation  of  Dr.  King’s  OMT  (King  1989;  King  and  Leon  1988).  In  OMT 
we  attempt  to  invert  the  equation  of  radiative  transfer 


R(P)=?W(p/P)B(p)^  (1) 

o 

where  R  is  the  observed  radiance  and  B  the  Planck  function.  It  is  convenient  to  parametrize 
this  problem  in  terms  of  the  atmospheric  pressure,  p ,  rather  than  altitude.  Equation  (1) 
represents  the  radiance  observed  in  a  single  channel  of  a  radiometer  measuring  the  upwelling 
radiance  at  the  top  of  the  atmosphere.  If'  is  a  weight  function  for  that  channel  and  P  is  a 
characteristic  pressure  value  for  that  channel,  crudely  speaking  the  centroid  of  the  weight 
function.  (Our  work  contains  significant  new  developments  in  the  definition  of  P ,  discussed  in 
Section  2.3.)  The  weight  function  is  defined  to  be 


W  « 


dr 

d  In  p 


(2) 


where  r  is  the  atmospheric  transmittance.  Thus  If'  measures  the  interaction  between  the 
radiation  field  and  atmospheric  material. 

It  is  often  convenient  to  parametrize  the  weight  function  in  terms  of  a  generalized 
exponential  weight  function 


Wm(p/P)  «  Wm(x)  =  i|~ jy  *  exp(-rnx1/m) 


(3) 


where  x  =  p/P  and  m  is  a  fitting  parameter  controlling  the  width  of  the  generalized  exponential 
weight  function,  and  T  denotes  the  usual  T  function  (King  1985;  King,  Hohlfeld,  and  Kilian 
1989).  This  simple  functional  form  captures  the  essential  physical  character  of  atmospheric 
weight  functions  in  many  cases.  (For  m  =  1  this  formula  reduces  to  a  simple  exponential 
weighting  function.)  While  OMT  bay  be  carried  out  with  completely  general  weight  functions, 
for  the  generalized  exponential  weight  functions  of  the  form  given  in  equation  (3)  the 
derivation  is  particularly  clear. 

When  we  substitute  the  generalized  exponential  weight  function  into  the  equation  of 
radiative  transfer,  we  obtain 


R(P)  =  nTSTTTp-C  B(p)  exPl‘m< P/P),/ml 


(4) 


For  m  -  I,  Eq.  (4)  immediately  assumes  the  character  of  a  Laplace  transform. 
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R(P)  =  L  j“B(p)e'(p/P) 


dp 


(5) 


i.e.  PR(P)  is  the  Laplace  transform  of  B(p)  with  respect  to  the  transform  variable  \/P.  When  m 
*  J,  the  radiance  profile  and  Planck  function  profile  are  then  related  by  a  generalization  of  the 
Laplace  transform  defined  in  terms  of  the  generalized  exponential  function.  All  of  the  useful 
analytic  properties  of  the  Laplace  transform  are  retained  in  this  generalization  of  the  Laplace 
transform  utilizing  the  generalized  exponential  function  kernel. 

If  a  choice  of  functional  form  is  made  for  B(p)  (or  alternatively  R(P)),  equation  (4)  then 
immediately  implies  the  functional  form  R(P)  (alternatively  B(p))  must  assume.  This 
observation  is  especially  pertinent  if  we  note  that  the  Planck  function  profile  of  a  radiative 
atmosphere  is  exponential  in  form,  i.e. 


B(p)  =  Le'kp 


(6) 


with  k  and  L  constants  (Chandrasekhar  1960).  The  corresponding  {m  =  1)  radiance  profile 
function  is  then. 


R(P)  -  fijjr  (7) 

This  discussion  motivates  a  choice  of  the  functional  form  by  which  radiance  data  is  represented 
(King  1989;  King  and  Leon  1988), 


R(P)  -  a  ♦  bP  +  £  ,-^p  (8) 

where  j  hyperbolic  terms  are  included.  The  addition  of  a  linear  term  of  form  a  +  bP  (with  a 
and  b  constants)  was  found  by  King  and  Leon  to  be  of  practical  utility  in  the  representation  of 
radiance  data.  Generalization  to  the  inclusion  of  a  polynomial  of  arbitrary  order  in  P  is 
straightforward.  The  corresponding  Planck  function  profile  obtained  by  generalized  inverse 
Laplace  transformation  (for  arbitrary  m)  of  Eq.  (7)  is 


B(p)  -  a  ♦  bp  +  t  LiEm(-kiP) 

»=i 


(9) 


Where  E Ax)  is  the  generalized  exponential  function  with  parameter  m.  This  new  function  is 
discussed  and  its  general  properties  exhibited  in  the  following  section. 

The  choice  of  an  exponential  (or  generalized  exponential  form  of  the  Planck  function 
profile,  motivated  by  the  expectation  of  the  real  atmospheric  acting  at  least  in  part  as  a 
radiative  atmosphere,  indicates  that  Eq.  (8)  is  a  natural  choice  of  functional  form  for  the 
representation  of  radiance  data. 
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Profiles.  In  this  section  we  briefly  describe  some  of  the  properties  of  the  generalized 
exponential  function  introduced  in  the  previous  section.  The  generalized  exponential  weighting 
function,  Em(x),  has  a  power  series  expansion 


^2 

£m(x)  =  (xj2 Y+  + 


+  ... 


(10) 


where  the  n  ’s  are  defined  in  terms  of  moments  of  the  generalized  exponential  weight  functions 
by 


*Un+D  =  Jfl  xn  Wm(x) 


dx 

m'- l'1  X~  ’ 


(ID 


It  is  apparent  that  equation  (10)  represents  a  generalization  of  the  usual  exponential  function, 
ex,  if  we  write  the  series  expansion  of  ex  using  F(n+1)  to  represent  «!, 


ex 


ncn+l)  + 


(12) 


and  note  that  T(«)  is  defined  by  its  integral  representation 


I\n)  =  t"’1  e  *  dt  ,  (13) 

and  so  £  (x)  =  ex.  This  motivates  the  form  of  equation  (9).  We  plot  in  Figure  1  a  family  of 
curves  of  Em(x )  for  various  values  of  m.  It  is  apparent  from  this  plot  that  variation  of  m 
controls  the  stiffness  of  the  resulting  Planck  function  curve. 

The  definition  of  the  generalized  exponential  function  discussed  here  is  motivated  by  the 
functional  form  of  the  generalized  weight  function  IFm(x).  However,  an  empirically  obtained 
weight  function  of  general  form  may  also  be  reduced  to  a  collection  of  moments  and  a 
corresponding  function  E(x)  computed,  and  so  this  development  of  OMT  is  not  restricted  simply 
to  the  use  of  generalized  exponential  weight  functions. 

Once  a  fit  has  been  made  to  the  radiance  data  on  the  basis  of  the  form  of  equation  (8), 
translation  into  the  form  of  the  Planck  function  profile  given  by  equation  (9)  is  immediate, 
provided  that  we  have  determined  m  values  for  the  evaluation  of  E  (-kp)  based  on  a  fit  to  the 
atmospheric  weight  functions  characterizing  the  instrumental  channels  used  in  the  fit.  A 
wavenumber  is  then  chosen  (usually  the  wavenumber  of  the  central  channel,  to  which  radiance 
measurements  are  normalized)  and  the  Planck  formula. 
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B 


(14) 


for  a  given  height  (parametrized  by  its  pressure  p),  may  be  solved  for  the  temperature  at  that 
height.  Here  v  is  the  wavenumber,  and  and  r,  are  constants  with  values,  =  1.19107  x  10'6 
ergcm2/secsr  and  c2  =  1.438858  cm-K,  and  T  in  degrees  Kelvin.  The  mathematical 
formalism  outlined  here  gives  a  clear  path  from  the  radiance  measurements  to  the  specification 
of  the  temperature  as  a  function  of  height  in  the  atmosphere.  It  is  important  to  note  that  the 
OMT  transformation  (of  equation  (8)  to  equation  (9))  and  the  solution  of  the  atmospheric 
temperature  profile  from  the  Planck  function  profile  are  one-to-one  transformations,  i.e.  given 
a  particular  choice  of  radiance  basis  functions,  as  in  equation  (8),  once  a  fit  is  made  to  that 
functional  form,  the  atmospheric  temperature  profile,  T  =  T(h ),  is  uniquely  determined. 


2.3  Definition  of  P  Values.  In  work  undertaken  on  Differential  Inversion  and  OMT  to  date 
(King,  Hohlfeld,  and  Kilian  1989)  P  has  been  defined  as  that  value  of  the  pressure  at  which  the 
weight  function  achieves  its  maximum  value.  We  have  noted  that  this  definition  is  not  entirely 
satisfactory  in  that  some  weight  functions  may  have  a  complicated  geometric  form,  and  may  in 
fact  exhibit  multiple  maxima,  particularly  if  several  absorbing  chemical  species  are  present. 

In  view  of  the  character  of  equation  (1),  the  equation  of  radiative  transfer,  which  specifies 
the  generation  of  the  upwelling  radiance  based  on  the  distribution  of  temperature  versus  height 
in  the  atmosphere,  we  have  begun  investigations  of  more  suitable  definitions  of  P  based  on 
moments  of  the  weight  function.  In  particular,  if  we  consider  equation  (1)  for  the  case  of  an 
isothermal  atmosphere,  B{p)  will  be  a  constant  function,  and  the  contributions  to  R(P)  will  be 
centered  about  the  center  of  mass  of  the  weight  function,  i.e.  the  centroid.  This  suggests  that 
the  centroid  of  the  weight  function  ought  to  be  chosen  as  a  value  for  P,  and  this  definition  is 
attractive  in  that  it  is  well  defined  even  for  weight  functions  which  have  complicated 
geometries  with  multiple  maxima.  Given  the  definition  of  the  weight  function,  as  in  equation 
(2),  it  is  apparent  that  this  definition  is  equivalent  to  selection  of  that  pressure  value  at  which  r 
=  1/2. 

At  present,  this  argument  is  suggestive,  rather  than  rigorous.  However,  we  have  noted  in 
our  research  efforts  attempting  to  construct  OMT  fits  to  TOVS  radiance  data,  that  superior  fits, 
in  the  sense  of  having  lower  x.2  values,  were  obtained  using  these  P  values  defined  in  terms  of 
the  centroid  of  the  TOVS  weight  functions.  Further  research  will  be  required  during  the  Phase 
II  effort  to  investigate  rigorously  appropriate  definitions  of  P.  This  problem  is  of  special 
significance  for  weight  functions  which  have  significant  values  at  ground  level  (and  are 
therefore  truncated  at  ground  level  with  significant  loss  of  the  area  under  the  weight  function 
curve). 


2.4  Development  of  the  Nonlinear  Least  Squares  Nonlinear  Hyperbolic  Algorithm  (NLLS 
NHA1 

2.4.1  Motivation.  Leon  and  King  (1988)  have  investigated  an  algorithm  for  directly  fitting  a 
formula  of  the  form  of  equation  (8)  to  radiance  data,  and  despite  the  attractive  adaptive 
features  of  that  algorithm,  the  requirements  on  the  accuracy  of  the  radiance  values  are  severe. 
This  property  arises  from  the  character  of  the  Leon  and  King  Nonlinear  Hyperbolic  Algorithm 
(NHA)  as  an  interpolation  formula;  the  solution  resulting  from  the  NHA  is  constrained  to  pass 
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identically  through  each  of  the  radiance  values.  In  practice,  radiance  data  are  corrupted  by 
noise  (usually  with  Gaussian  statistics)  and  the  requirement  of  passing  through  each  data  point  is 
unnecessarily  stringent  under  some  circumstances. 


Furthermore,  fitting  a  formula  of  the  form  of  equation  (8)  is  a  difficult  numerical  problem. 
Fitting  a  sum  of  exponential  functions  to  a  collection  of  data  is  a  classical  example  of  an  ill- 
conditioned  fitting  problem  (Acton  1970).  This  problem  is  very  similar  because  e"x  and  1/(1  + 
x )  have  very  similar  shapes  over  the  first  octave  of  their  range.  Consequently,  this  problem 
may  be  expected  to  exhibit  similar  ill-conditionedness,  as  is  in  fact  observed. 


Motivated  by  these  considerations,  we  have  begun  development  of  a  new  algorithm  for 
fitting  radiance  data  to  a  formula  of  the  form  of  equation  (1)  based  on  nonlinear  least  squares 
techniques  (Hohlfeld,  Kilian,  Drueding,  and  Ebersole  1988).  We  have  named  this  algorithm  the 
NLLS  NHA.  The  intent  of  the  NLLS  NHA  is  to  fit  as  many  terms  of  the  formula  of  equation 
(1)  as  can  be  supported  by  the  radiance  data  set,  as  determined  by  objective  criteria  of  the 
“goodness  of  fit",  such  as  measures  of  x2-  Ultimately,  based  on  the  development  of  NLLS 
NHA  techniques,  it  will  be  possible  to  derive  measures  of  the  information  content  of  a  radiance 
data  set  based  on  the  number  of  well  conditioned  terms  of  the  fit  which  can  be  obtained. 

Nonlinear  least  squares  fitting  algorithms  operate  by  the  minimization  of  an  appropriately 
constructed  cost  function,  usually  termed  x2 ■  defined  by 


X 


2 


£  f  y.-><a1,...ak,xi) 

L  a 

i=i  L. 


(14) 


The  algorithm  seeks  to  minimize  x2.  by  varying  the  k  parameters  of  the  theory.  a . ak ,  and  a 

total  of  n  measurements  have  been  taken  at  n  values  of  independent  variable,  x.,  to  obtain 
values  of  the  dependent  variable,  y .  The  description  of  the  model  is  contained  in  the  function 
y(av...ak,x). 


Choice  of  the  appropriate  algorithm  for  achieving  x2  minimization  in  a  given  functional 
fitting  problem  is  most  readily  made  on  the  basis  of  an  investigation  of  the  geometric  properties 
of  the  x2  surface.  The  information  required  for  making  such  a  choice  of  algorithm  is  described 
below. 

2.4.2.  Selection  of  the  Levenbere-Marauardt  Algorithm  for  the  NLLS  NHA.  Work  undertaken 
earlier  in  this  project  has  characterized  the  geometry  of  the  x2  surface  generated  by  attempting 
to  fit  a  functional  form  of  the  radiance  data,  as  given  by  equation  (8)  to  radiance  data  from  the 
TIROS  Operational  Vertical  Sounder  (TOVS).  We  have  found  that  the  characteristic  size  of  the 
minima  in  a  and  the  L' s  for  the  positive  k  portion  of  the  radiance  x2  surface  is  typically  of  the 
order  of  a  few  tens  in  magnitude  and  the  minima  in  the  A/s  are  of  the  order  of  unity  in 
characteristic  width.  Minima  in  the  radiance  x2  surface  with  Ac  <  0  (see  section  3.1)  are  much 
more  constrained  in  L  (width  <  10)  and  in  k  (width  <  1/2).  Multiple  deep  minima  in  x2  exist  in 
close  proximity  in  this  region.  (Sec  Figure  2.)  The  geometry  of  the  x2  surfaces  in  the 
immediate  neighborhood  of  the  minima  is  regular  and  the  ellipses  defined  by  levels  of  constant 
X2  are  reasonably  well  aligned  from  one  level  to  the  next. 

In  Figure  3  we  show  the  x2  surface  generated  by  comparing  the  temperature  profile 
corresponding  to  the  radiances  generating  Figure  2  against  the  corresponding  radiosonde 
observations.  It  is  important  to  note  that  the  minimum  in  the  temperature  x2  surface  coincides 
with  the  global  minimum  in  the  radiance  x2  surface  shown  in  Figure  2,  thus  confirming  our 
identification  of  the  correct  physical  solution  in  the  OMT  processing  of  this  radiance  scan. 
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Parameter  3 


ONTO'JR  VALUES 


2  *  +4.040E  +  01 

3  *  +1  .927E  +  02 

4  *  +  9.200E  +  02 

6  *  +2.095E  +  04 


Figure  2:  y2  surface  of  the  Optical  Measure  Theory  fit  to  a  radiance  profile.  Shown  here  is  a 
contour  map  of  the  x2  surface  generated  using  OMT  with  four  parameters  to  fit  a  radiance 
profile.  The  TOVS  radiance  scan  from  Zone  1  Pass  9  was  chosen  for  this  example.  Two 
parameters  are  held  constant  (  a  ■  20  and  b  =  45)  and  two  are  varied  (  Ki  ranging  from  -5  to  5 
and  Li  ranging  from  -50  to  50).  The  geometry  seen  here  is  representative  of  all  the  TOVS 
scans. 
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Figure  3:  Surface  of.  the.  Optical  Measure  Theory  fit  to  the  temperature  profile  from  RAOB 
data.  Shown  here  is  the  contour  map  of  the  x3  surface  generated  using  OMT  to  fit  a 
temperature  profile  from  RAOB  data  under  the  same  conditions  as  Figure  2. 
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Owing  to  this  comparatively  well-behaved  geometry  in  the  x2  surfaces,  we  are  motivated  to 
select  a  comparatively  conventional  but  well-tested  algorithm  for  the  x2  minimization  in  the 
NLLS  NHA,  the  Levenberg-Marquardt  algorithm.  This  algorithm  is  discussed  by  Press, 
Flannery,  Teukolsky,  and  Vetterling  (1986),  and  we  reproduce  some  relevant  elements  of  their 
discussion  in  the  following  section  to  motivate  our  choice  of  this  algorithm  for  the  NLLS  NHA. 

2.4,3  Development  of  the  Levenbere-Nlarauardt  Algorithm.  We  consider  here  the  general 

case  when  a  model  depends  nonlinearly  on  a  set  of  M  unknown  parameters  ak,  k  *  1,2 . M. 

[The  parameters  a ^  correspond  to  the  parameters  a,  6,  k},  Ly  k2  and  L%  in  equation  (9).]  Our 
objective  is  to  minimize  the  x2  merit  function,  as  defined  in  equation  (14),  and  thereby  to 
determine  the  best-fit  parameters  for  the  fit  of  our  model.  Since  the  model  is  nonlinear,  it  is 
necessary  to  carry  out  an  iterative  solution  for  our  desired  parameters.  Given  trial  values  for 
the  parameters,  some  procedure  improving  the  trial  solution  is  carried  out  until  x2  stops  (or 
effectively  stops)  decreasing. 

Sufficiently  close  to  the  minimum,  we  can  expect  that  the  x2  function  will  be  approximated 
by  a  quadratic  form,  which  may  be  written  as 


X2(a)  =  7  -  d  a  +  J[-a>D-a  (16) 

where  d  is  an  M-vector  and  D  is  an  M  x  M  matrix.  If  the  approximation  is  a  good  one,  we  can 
jump  from  the  current  trial  parameters  acur  to  the  minimizing  ones  amin  in  a  single  step, 


=  a. 


*  »  1  [-vx8(>,„,)  ] 


07) 


On  the  other  hand,  equation  (17)  might  be  a  poor  local  minimization  to  the  shape  of  the 
function  we  are  attempting  to  minimize  at  acur.  In  such  an  instance,  all  that  can  be  done  is  to 
take  a  step  along  the  gradient  of  x2,  as  is  done  in  the  steepest  descent  method  for  minimization. 
In  that  case  the  next  iteration  would  be  given  by 


anext  =  acur  -  constant  x  Vx2(acur) 


(18) 


where  the  constant  is  chosen  small  enough  so  as  to  not  exhaust  the  downhill  direction. 

To  use  either  equation  (17)  or  (18),  we  must  be  able  to  compute  the  gradient  of  the  x2 
function  at  any  set  of  parameters  a.  To  use  equation  (4),  we  also  need  the  matrix  D,  the  second 
derivative  matrix  (Hessian  matrix)  of  the  x2  merit  function,  at  any  a.  In  the  Levenberg- 
Marquardt  method,  we  shall  use  either  the  Hessian  matrix  or  the  method  of  steepest  descents  for 
the  minimization  of  x2»  depending  upon  the  information  that  can  be  extracted  by  the  algorithm 
concerning  the  local  geometry  of  the  x2  surface. 
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The  model  to  be  fitted  is 


y  =  y(x;a) 


and  the  v2  merit  function  is  given  by 


,  n  r y i-yi***)  l2 

*  <*>  ■  s  L  "i  J 


The  gradient  of  x2  with  respect  to  the  parameters  a,  which  will  be  zero  at  the  x2  minimum,  has 
components 


gx2  ^  n  [y.-yfx^j  ^(x.;a) 

*»k  "  k  °i 2  ^k 

Taking  an  additional  partial  derivative  yields 


f  L  Iv  .  ^ VA.  1 

aak0a,  ~  k  of  L  ^k  *1  Iyi  y(xi’a;j  aa,aak  J 


It  is  conventional  to  remove  the  factors  of  2  by  defining 


a  _  I  aA 
~  fin 


Q  =  L  fxL 

w  "  2  3au5a. 


making  [aj  =  (1/2)D  in  equation  (17),  in  terms  of  which  that  equation  can  be  rewritten  as  the 
set  of  linear  equations 


I  Qkl  *|  *  4 
1  =  1 
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This  set  of  equations  is  solved  for  the  increments  5a,,  which  are  added  to  the  current 
approximation  to  obtain  the  next  approximation  to  the  set  of  parameters.  In  the  context  of 
least-squares  fitting,  the  matrix  [a],  which  is  equal  to  one-half  times  the  Hessian  matrix,  is 
usually  called  the  curvature  matrix. 

Equation  (18),  the  steepest  descent  formula  may  be  expressed 


5a,  =  constant  x  f) 


(26) 


Note  that  the  components  ak,  of  the  Hessian  matrix  depend  both  on  the  first  derivatives  and 
on  the  second  derivatives  of  the  basis  functions  with  respect  to  their  parameters.  Many 
treatments  proceed  to  ignore  the  second  derivative,  without  justification,  which  we  shall  attempt 
to  supply  here.  Second  derivatives  occur  because  the  gradient  [equation  (21)]  already  has  a 
dependence  of  dy./da.,  and  so  the  next  derivative  must  simply  contain  terms  involving 
d^y.Jda^da, .  The  second  derivative  term  can  be  neglected  when  it  is  zero  (as  in  linear 
regression),  or  small  enough  to  be  negligible  when  compared  to  the  term  involving  the  first 
derivative.  It  also  has  an  additional  possibility  of  being  ignorably  small  in  practice:  The  term 
multiplying  the  second  derivative  in  equation  (22)  is  [y,  -  >-(x.;a)J.  For  a  successful  model,  this 
term  should  be  just  the  random  measurement  error  of  each  point.  This  error  can  have  either 
sign,  and  should  in  general  be  uncorrelated  with  the  model.  Therefore,  the  second  derivative 
terms  tend  to  cancel  out  when  summed  over  /. 

Inclusion  of  the  second-derivative  term  can  in  fact  be  destabilizing  if  the  model  fits  badly 
or  is  contaminated  by  outlier  points  that  are  unlikely  to  be  offset  by  compensating  points  of 
opposite  sign.  From  this  point  on  we  shall  always  use  as  the  definition  of  ak,  the  formula 


£  i  ray(x>;a)  1 


(27) 


It  should  be  understood  that  minor  (or  even  major)  fiddling  with  [a]  has  no  effect  at  all  on 
which  final  set  of  parameters  a  is  reached,  but  only  affects  the  iterative  route  that  is  taken 
getting  there.  The  condition  at  the  x2  minimum,  that  0k  =  0  for  all  k,  is  independent  of  how 
[a]  is  defined. 

2.4.4  The  Levenbere-Marauardt  Method.  An  elegant  method  has  been  put  forward  by 
Marquardt  (1963)  based  on  an  earlier  suggestion  bv  Levenberg,  which  combines  smoothly  the 
inverse-Hessian  and  steepest  descent  methods  for  x2  minimization,  as  discussed  in  the  previous 
section.  The  steepest  descent  method  is  used  far  from  the  x2  minimum,  when  the  geometry  of 
the  x2  *s  not  well  modeled  by  a  locally  quadratic  behavior.  Then,  as  the  minimum  is 
approached,  the  Levenberg-Marquardt  algorithm  is  switched  continuously  to  x2  minimization  by 
the  inverse-Hessian  method.  The  Levenberg-Marquardt  method  works  well  in  practice,  and  has 
become  a  standard  nonlinear  least-squares  algorithm  (Press,  Flannery,  Teukolsky,  and  Vetterling 
1986). 

Application  of  the  Levenberg-Marquardt  method  depends  upon  two  observations  regarding 
the  steepest  descent  and  inverse-Hessian  methods.  Consider  the  constant  multiplier  in  equation 
(26);  within  the  context  of  a  standard  steepest  descent  implementation,  we  have  no  information 
regarding  the  magnitude  of  the  constant  or  about  the  scale  over  which  it  may  be  applied.  No 
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information  answering  these  questions  is  available  in  the  calculation  of  the  gradient  of  x2>  which 
tells  us  only  the  value  of  the  slope  of  the  x2  surface,  not  the  distance  that  the  slope  extends. 
Marquardt’s  first  insight  is  that  the  components  of  the  Hessian  matrix,  even  if  they  are  not 
usable  in  any  precise  fashion,  give  some  information  regaiding  the  order-of-magnitude  scale  of 
the  problem. 

The  quantity  x2  is  dimensionless,  i.e.  is  a  pure  number.  On  the  other  hand,  0k  has  the 
dimensions  of  l/ak,  which  may  be  a  dimensional  quantity  (and  in  fact  each  component  of  0k 
may  have  different  dimensions).  The  constant  of  proportionality  between  £k  and  6a.  must 
therefore  have  the  dimensions  of  a  2  In  the  components  of  [a],  the  obvious  quantity  with  these 
dimensions  is  l/akk,  the  reciprocal  of  the  diagonal  element,  which  must  be  the  scale  of  the 
constant  for  the  application  of  steepest  descent.  In  practice,  it  is  necessary  to  divide  the 
constant  by  some  (nondimensional)  fudge  factor  A,  with  the  possibility  of  setting  A  »  1  to  cut 
down  the  step  size.  In  other  words,  replace  equation  (26)  with 


-  4  f,  <28> 

It  is  necessary  that  be  positive,  but  this  is  guaranteed  by  our  definition  of  equation  (27), 
which  was  one  of  the  motivations  for  making  that  choice. 

The  second  insight  in  Marquardt’s  development  of  this  algorithm  is  that  equations  (28)  and 
(28)  can  be  combined  if  we  define  a  new  matrix  q  by  the  following  prescription 


\  s  Qjj  < 1  +  A) 

7jk  s  Qjit  (j  *  k)  (29) 

and  then  replace  both  equations  (25)  and  (28)  by 

£  =  K  (30) 

1=1  Kl  * 


When  A  is  very  large,  the  matrix  -j  is  forced  into  being  diagonally  dominant,  so  that  equation 
(30)  goes  over  to  the  form  of  equation  (28).  On  the  other  hand,  as  A  approaches  zero,  equation 
(30)  goes  over  to  the  form  of  equation  (25).  What  has  been  accomplished  by  use  of  j  in  the  x2 
Levenberg-Marquardt  algorithm  is  to  have  the  minimization  algorithm  convert  smoothly  from 
the  steepest  descent  algorithm  to  the  inverse-Hessian  minimization  algorithm  as  the  minimum  in 
the  x2  surface  is  approached. 


2.4.5  Limitations  of  the  Levenbere-Marauardt  Method.  Our  development  of  the  NLLS  NHA 
has  reduced  the  process  of  obtaining  OMT  solutions  to  the  process  of  minimization  on  an 
appropriately  defined  x2  surface,  which  we  accomplish  by  the  Levenberg-Marquardt  method. 
However,  while  the  Levenberg-Marquardt  method  is  a  relatively  sophisticated  minimization 
technique,  it  still  takes  account  only  of  locally  defined  information  on  the  x2  surface,  and  each 
step  in  parameter  space  is  made  so  as  to  reduce  the  value  of  x2  of  the  solution.  Accordingly,  it 
is  possible  for  the  solution  to  become  “trapped”  in  a  local  minimum  in  the  x2  surface,  and  thus 
never  attain  the  true  global  solution.  We  circumvent  this  problem  in  our  present  research  by 
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selection  of  appropriate  starting  vectors,  motivated  by  our  analysis  of  the  surface,  which 
permit  us  to  avoid  the  local  minima,  while  stilt  not  biassing  the  solution.  This  procedure  works 
with  reasonably  high  probability,  but  is  still  somewhat  ad  hoc,  and  a  more  theoretically 
satisfactory  alternative  should  be  sought. 

Simulated  annealing  (Metropolis  cl  «/.  1953)  and  neural  network  techniques  (Hopfieid  and 
Tank  1985)  have  been  applied  to  minimization  problems  in  radiative  transfer  (Jeffrey  and 
Rosner  1986).  These  techniques  have  been  shown  to  yield  true  global  minima  with  effectively 
unit  probability.  While  these  techniques  are  computationally  intensive,  being  implemented  most 
effectively  in  parallel  computer  architectures,  they  are  of  great  interest  and  have  potentially 
significant  applications  in  the  inversion  problem. 
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3.  ANALYSIS 


3.1  New  QMT  Solutions  Our  previous  OMT  solutions  (Hohlfeld,  Kilian,  Drueding,  and 
Ebersole  1988)  exhibited  the  character  of  being  a  reasonable  fit  at  the  tropopause  and  for  the 
temperature  profile  in  the  stratosphere  (in  the  sense  of  fitting  global  morphological  features  of 
the  temperature  profile,  rather  than  accounting  for  all  the  details  of  the  radiosonde 
observations).  However,  these  initial  OMT  fits  deviated  significantly  from  radiosonde 
measurements  of  the  atmospheric  temperature  in  the  troposphere,  by  an  amount  increasing  as 
the  pressure  increases,  up  to  some  tens  of  Kelvins  at  ground  level.  All  such  fits  contained  a 
single  hyperbolic  term,  with  a  positive  value  for  A 

Subsequent  analysis  has  shown  that  a  negative  value  for  Aj  yields  significantly  improved  fits 
to  both  the  radiance  profile  and  temperature  profile  (as  measured  by  the  x2  calculation  in  the 
NLLS  NHA),  particuiarly  over  the  domain  of  tropospheric  pressures.  A  value  of  Aj  less  than 
zero  may  initially  seem  troubling  theoretically  in  that  it  corresponds  to  a  pole  in  the  OMT 
radiance  formula 


R(P)  -  a  +  bP  ♦  £  ,-J^p  (31) 

at  positive  P  ( i.e .  at  positive  pressure).  Here  P  is  the  pressure  of  the  maxima  of  the  radiance 
channels,  suitably  normalized  (to  a  pressure  of  1  bar  for  this  work,  so  that  0  <  P  <  1  over  the 
physical  atmosphere).  Values  of  A\  such  that  -1  <  A.  <  0,  correspond  to  poles  in  R(n)  occurring 
below  ground  level,  i.e.  in  an  unobservable  region 

With  the  viewpoint  that  equation  (1)  constitutes  a  rational  function  representation  of  the 
radiances  over  the  observed  pressure  range,  (he  existence  of  a  pole  in  the  radiance  formula  for  a 
value  of  the  pressure  at  which  radiance  measurements  cannot  be  obtained  is  not  a  cause  for 

concern. 

We  note  that  in  the  present  formulation  of  OMT,  some  difficulties  arise  due  to  the 
truncation  of  weight  functions  at  ground  level  (P  =  1).  This  prompts  the  speculation  (albeit  a 
reasonable  one)  that  the  negative  A:  values  obtained  in  our  OMT  solutions  arise  due  to  this 
truncation  of  the  atmospheric  weight  functions.  This  motivates  an  extension  of  OMT,  beyond 
its  present  formulation,  which  more  realistically  treats  the  radiative  transfer  problem  at  the  base 
of  the  atmosphere.  Further  research,  undertaken  during  the  Phase  II  effort,  should  clarify  these 
important  practical  and  theoretical  issues. 


3.2  Description  of  NLLS  NHA  Software. 

3.2. J  Introduction.  We  have  designed  a  computerized  algorithmic  framework  to  test  different 
implementations  of  x2  minimization  algorithms.  This  framework  provides  several  options  when 
running  a  algorithm.  We  look  for  a  set  of  options  which  produce  a  good  fit  to  the  radiance  data 
and  are  applicable  to  all  available  cases.  The  algorithms  for  the  Non-Linear  Least  Squares,  the 
Non-Linear  Hyperbolic  fitting,  and  a  x2  analysis  of  a  fit  to  radiance  data  are  written  in 
FORTRAN.  We  created  a  shell  program  to  drive  these  routines,  to  report  the  results,  and  to 
quickly  change  the  methods  of  implementation.  The  program  runs  in  two  modes,  a  highly 
interactive  mode  where  we  chose  one  scan  and  make  immediate  adjustments  to  our  method  and 
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a  batch  mode  which  analyzes  all  scans  using  one  met'  1.  The  batch  mode  also  provides  the 
means  of  processing  the  data  when  we  determine  an  optima!  robust  reduction  procedure  using 
these  algorithms.  A  spread  sheet  program  (using  Microsoft  Excel1*0)  gives  a  more  detailed 
analysis  of  the  radiance  fits,  and  presides  graphs  of  such  fits.  The  spread  sheet  program  does 
not  run  the  NHA  or  N'LLS  algorithms;  it  requires  information  on  tne  constants  in  OMT  used  to 
represent  the  radiance  curve.  Using  the  Microsoft  Windows®  environment  we  run  both  the 
i  OR  I  RAN  program  Un  interactive  mode),  and  tne  spread  sheet  program  simultaneously  to 
immediately  analyze  in  detail  the  results  of  running  Nl.LS-NHA  algorithms  and  to  view  the 
graphs  of  the  ot  the  OMT  lits  as  the  adjusinuufs  a;-,  made  to  our  methods. 


3.2.2  FORTRAN  Program  We  chose  the  Levenberg-Marquardt  method  for  the  NLLS 
algorithm  (Press,  Flannery,  Teukolsky  and  Vetterling,  1 9S6).  We  took  the  basic  algorithm  from 
Press  et  al.  and  integrated  it  with  the  shell  program,  BP  IT.  BUT  breaks  into  six  independent 
parts;  (1)  radiance  data  input,  (2 »  initialize  MIA  fitting,  (3)  run  NHA  fitting,  (4)  initialize 
NLl  S  fitting,  (5)  run  NLLS  fitting,  and  (6)  output  solution  parameters  and  fitted  radiance 
profiles. 

The  input  routine  reads  a  data  file  winch  contains  information  on  a  particular  TOVS  scan. 
Each  TOVS  scan  provides  seven  pairs  of  P  oar  and  Radiance  values.  When  BFIT  is  run 
interactively,  we  choose  one  TOVS  scan  to  test  in  batch  mode  BFIT  loops  though  all  scans. 

We  integrated  the  original  Non -I. meat  Hyperbolic  fitting  Algorithm  to  finish  thoroughly 
analyzing  that  method  and  its  usefulness.  When  tunning  t fie  NHA  routine,  we  choose  the 
number  of  constants  to  fit  te.g  ,  a  choice  of  four  constants  determines  one  constant  term,  one 
linear  term,  and  one  hyperbolic  term)  and  then  decide  which  channels  to  fit.  The  TOVS  data 
has  a  total  of  seven  channels. 

Given  a  number  of  channels,  (he  NHA  fitting  Algoiithm  provides  a  exact  fit  to  those 
channels.  1  he  algorithm  is  numerically  unstable  owing  to  its  character  as  an  interpolation 
algorithm.  Consequently,  it  piovides  very  different  solutions  depending  on  which  channels  are 
chosen,  and  it  cannot  determine  a  solution  for  more  than  four  channels  of  TOVS  data.  The 
NLLS  NHA  avoids  these  difficulties  because  it  is  not  an  interpolation  algorithm. 

BFIT  runs  the  NLLS  routine  according  to  3  script.  1  he  script  determines  which  constants 
are  to  be  varied  and  in  "  hat  order.  For  exam,  !e  we  choose  a  starting  vector  close  to  zero 
(presumably  a  bias-free  choice)  and  we  decide  to  release  the  constant  and  linear  term  only.  The 
then  algorithm  fits  the  data  with  no  hyperbolic  terms.  Then  we  release  one  hyperbolic  term  to 
see  how  it  corrects  the  fit.  The  algorithm’s  performance  also  depends  on  a  number  of 
parameters  such  as  the  estimated  errors  in  the  data  which  weight  different  channels  and  the 
derivatives  of  the  OMT  constants  which  weight  different  terms  in  the  curve,  BFIT  sets  these 
parameters. 

The  NLLS  is  run  according  to  the  script,  and  it  attempts  to  reduce  the  global  x2  fit  to  the 
radiance  data  by  varying  the  released  constants.  It  provides  what  it  determines  as  the  best  set 
o<  OMT  constants  for  the  fit. 

BFIT  reports  the  final  OMT  cm  .dams  and  provides  a  x2  Per  degree  of  freedom  analysis  of 
the  OMT  radiance  curve  fit  to  the  or  dinal  radiance  data.  This  information  is  logged  in  a  data 
file  to  maintain  a  record  of  each  method'  olution  and  performance. 
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3.2.3  Spread  Sheet  Program.  The  spread  sheet  program  analyzes  one  scan  in  detail.  All 
information  on  the  scan  is  extracted  from  the  data  files  used  by  the  FORTRAN  routines  and 
presented  in  clear  manner  in  one  work  sheet.  The  sheet  contains  the  location  and  time  of  the 
scan,  and  lists  in  table  form  all  information  about  each  channel  and  the  OMT  curve  match  to 
that  channel.  The  spread  sheet  program  does  several  different  error  calculations  in  the  tables. 

We  manually  enter  values  for  the  OMT  constant  used  to  create  the  radiance  curve.  If  we 
change  the  OMT  constants  on  the  work  sheet,  the  spread  sheet  program  immediately  recalculates 
the  OMT  radiance  curve  and  the  error.  A  graphics  sheet  takes  information  from  the  work  sheet 
and  creates  a  graph  of  the  radiance  data  points  and  the  OMT  curve.  The  graph  is  automatically 
updated  when  any  changes  are  made  on  the  work  sheet.  The  spread  sheet  program  quickly 
presents  clear  and  detailed  information  on  the  OMT  fit  but  does  not  provide  information  on 
how  the  fit  was  determined. 


3.3  Reduction  of  TOYS  Data  to  Atmospheric  Temperature  Profiles  using  OMT.  We  have 
readily  available  a  set  of  45  TOVS  radiance  scans,  15  from  each  of  the  three  latitude  zones. 
Each  scan  has  corresponding  radiosonde  observations  for  comparison  with  the  OMT  temperature 
profile.  These  scans  have  been  processed  with  the  NLLS  NHA  software  discussed  in  previous 
reports.  (See,  for  example,  Hohlfeld,  Kalian,  Drueding,  and  Ebersole  1988.) 

Of  the  total  of  45  TOVS  scans  processed,  the  complete  set  has  yielded  physically  meaningful 
solutions  involving  a  single  hyperbolic  term  (with  coefficients  ky  and  L.)  and  a  linear  term  with 
coefficients  a  and  b.  Results  of  these  runs  are  included  in  the  Appendix.  In  each  of  the  cases 
shown,  values  of  k{  were  obtained,  satisfying  -  I  <  k{  <  0,  which  corresponds  to  the  new  class  of 
solutions  discussed  in  the  previous  section. 

Low  latitude  atmospheric  temperature  profiles  exhibit  a  sharply  defined  tropopause  feature, 
as  evidenced  by  the  radiosonde  observations  for  the  TOVS  scans  which  we  have  studied.  It  is 
unrealistic  to  expect  OMT  fits  based  on  only  4  terms  to  adequately  represent  such  rapid 
variations  of  temperature  with  altitude.  We  plan  to  investigate  whether  this  situation  may  be 
partially  ameliorated  by  addition  of  a  second  hyperbolic  term,  or  the  addition  of  some  other 
term  to  equation  (8). 

We  have  obtained  meaningful  fits  to  both  long  wavelength  (667  cm1)  and  short  wavelength 
(2250  cm1)  TOVS  data.  However,  as  TOVS  has  only  5  short  wavelength  channels,  we  regard 
these  bits  as  being  less  significant  statistically  than  the  long  wavelength  OMT  fits. 

Examples  of  NHA  NLLS  processing  of  TOVS  data  to  OMT  radiance  and  temperature 
profiles  are  given  in  Figures  4  through  7.  For  the  TOVS  data  we  have  found  that  m  «  I  is 
appropriate  for  representing  TOVS  weight  functions  and  the  corresponding  representations  of 

3.4.  Discussion  of  the  Results  of  TOVS  Data  Reduction.  We  have  demonstrated  that  the  NLLS 
NHA  generates  OMT  fits  on  100%  of  a  set  of  45  TOVS  radiance  scans  from  all  latitude  zones. 
Thus  NLLS  NHA  is  a  robust  algorithm  suitable  for  processing  TOVS  data  on  a  routine  basis. 
Further  advances  in  the  construction  of  OMT  fits,  viewed  as  minimization  problems  in 
multidimensional  parameter  spaces,  will  allow  more  effective  and  efficient  bias-free 
computation  of  OMT  solutions. 
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Figure  4:  Optical  Measure  Theory  radiance  fit  to  TOYS  data  and  temperatuie  fit  to  RAOB  data. 
This  is  a  representative  case  of  the  fits  achieved  in  Zone  l  by  the  Non-Linear  Least  Squares 
algorithm  using  OMT  with  four  parameters  for  long  wavelength  data.  The  Non-linear  least 
squares  algorithm  was  given  an  initial  vector  with  a  =  b  »  Ki  »  0  and  Li  -  10  and  gave  a 
solution  with  a  *  20.23,  b  -  45.93,  Ki  -  0.07889  and  Li  *  12.29  for  Zone  I  Pass  9.  This  scan 
in  Zone  1  was  chosen  because  RAOB  data  was  available  up  to  10  mb. 


This  is  a  representative  case  of  the  fits  achieved  in  Zone  2  by  the  Non-Linear  Least  Squares 
algorithm  using  OMT  with  four  parameters  for  long  wavelength  data.  The  Non-linear  least 
squares  algorithm  was  given  an  initial  vector  with  a  =  b  »  Ki  *  0  and  Li  «  10,  and  gave  a 
solution  with  a  =  19.98,  b  =  85.54,  Ki  =  0.1233  and  Lt  ■  16.29  for  Zone  2  Pass  3.  This  scan  in 
Zone  2  was  chosen  because  RAOB  data  were  available  up  to  10  mb. 
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Figure  6:  Optical  Measure  Theory  radiance  fit  to  TOYS  data  and  temperature  fit  to  RAOB  data. 
This  is  a  representative  case  of  the  fits  achieved  in  Zone  3  by  the  Non-Linear  Least  Squares 
algorithm  using  OMT  with  four  parameters  for  long  wavelength  data.  The  Non-linear  least 
squares  algorithm  was  given  an  initial  vector  with  a  -  b  -  Ki  -  0  and  Li  ■  10  and  gave  a 
solution  with  a  -  15.40,  b  -  1 17.7,  Ki  -  0.008709  and  Li  -  18.49  for  Zone  3  Pass  12.  This  scan 
in  Zone  3  was  chosen  because  RAOB  data  were  available  up  to  10  mb. 
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TOVS  Radiance  Plot  (Short  wave  length) 
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Figure  7:  Optical  Measure  Theory  radiance  fit  to  TOVS  data  and  temperature  fit  to  RAOB  data. 
This  is  a  representative  case  of  the  fits  achieved  in  Zone  1  by  the  Non-Linear  Least  Squares 
algorithm  using  OMT  with  four  parameters  for  long  wavelength  data.  The  Non-linear  least 
squares  algorithm  was  given  an  initial  vector  with  a  =  b  -  Ki  *  0  and  Li  -  0.1  and  gave  a 
solution  with  a  »  -0.02927,  b  ■*  0.4591,  Ki  ■  1.1945  and  Li  ■  0.07186  for  Zone  1  Pass  9.  This 
scan  in  Zone  1  was  chosen  because  RAOB  data  were  available  up  to  10  mb. 
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It  is  highly  interesting  and  scientifically  significant  that  a  meaningful  OMT  fit  to  the 
radiance  data  which  yields  a  meteorological  reasonable  temperature  profile  can  be  constructed 
with  just  4  parameters.  First  of  all,  this  argues  that  the  choice  of  functional  form  in  equation 
(8)  was  an  appropriate  choice  of  the  representation  of  radiance  data.  Secondly,  it  is  highly 
suggestive  that  the  actual  information  content  of  the  TOVS  radiance  data  is  approximately  equal 
to  that  which  would  be  generated  by  4  statistically  independent  channels  (with  relative  errors 
characteristic  of  the  TOVS  radiometer).  Such  considerations  are  highly  significant  for  the 
design  of  atmospheric  sounders. 

We  note  that  the  atmospheric  temperature  profiles  obtained  using  the  NLLS  NHA  exhibit  a 
meteorological  character  and  have  a  smooth  form  following  the  general  trend  of  the  radiosonde 
observations.  The  solutions  we  have  obtained  with  -1  <  k.  <  0  are  a  significantly  better  match 
to  the  tropospheric  temperatures  than  our  previous  set  or"  solutions  obtained  for  k.  >  0.  It 
cannot  be  expected  that  an  OMT  temperature  profile  constructed  from  a  4  parameter  fit  to  the 
radiance  data  can  match  the  detailed  structure  of  the  radiosonde  observations.  (In  any  event, 
the  radiosonde  observations  available  to  us  from  NOAA  may  be  as  much  as  3  degrees  in  latitude 
and  longitude  from  the  geographic  location  of  the  TOVS  radiance  observations.)  However,  it 
can  be  seen  that  we  have  captured  the  overall  character  of  the  actual  atmospheric  temperature 
profile  in  the  OMT  temperature  profile.  Further  research  will  allow  interpretation  of  the  OMT 
fit  parameters  in  more  direct  meteorological  terms,  such  as  tropopause  temperature  and  pressure, 
and  temperature  lapse  rates. 
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4.  DIRECTIONS  FOR  FUTURE  RESEARCH 

The  research  undertaken  during  the  Interim  Research  Period  has  demonstrated  the  utility  of 
viewing  the  construction  of  OMT  fits  as  minimization  problems  in  a  multidimensional  parameter 
space,  and  has  shown  that  these  fits  exhibit  the  meteorological  characteristics  desired. 
Furthermore,  these  algorithms  are  robust,  working  on  our  entire  set  of  available  TOVS  radiance 
data.  This  represents  an  important  milestone  in  the  process  of  developing  operational  OMT 
capabilities.  Further  research  tasks  directed  toward  achieving  that  long-range  objective  can  now 
be  listed. 

o  Development  of  algorithms  for  global  minimization  on  the  radiance  x2  surface  utilizing 
nonlocal  information  on  x2  geometry  and  allowing  construction  of  nonlinear  cost  functions. 

•  Elucidation  of  the  mathematical  relationship  between  the  parameters  of  the  OMT  fit  and 
atmospheric  stiuctuie  parameters. 

•  Investigation  of  new  data  sets,  so  as  to  establish  that  the  results  we  have  obtained  are  not 
unique  to  the  TOVS  data  set.  Investigation  of  data  sets  with  larger  numbers  of  channels  is 
particularly  crucial. 

•  The  theoretical  formulation  of  OMT  should  be  extended  to  include  properly  the  effects 
of  the  ground,  both  in  terms  of  its  radiance  and  in  terms  of  the  truncation  of  atmospheric 
weight  functions. 

•  Analysis  of  OMT  algorithms  should  be  undertaken  from  information  theoretic  principles 
to  determine  the  information  content  of  available  radiance  data  sets  and  its  bearing  of  the 
determination  of  the  vertical  temperature  structure  of  the  atmosphere. 

In  conjunction  with  the  rest  of  our  proposed  Phase  II  research  effort,  considerable  progress 
can  be  made  leading  to  operational  implementation  of  OMT  algorithms. 
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We  have  demonstrated  an  implementation  of  Dr.  King's  Optical  Measure  Theory  (OMT) 
based  on  mathematical  concepts  of  x2  minimization  on  multidimensional  parameter  spaces.  This 
approach  is  embodied  in  our  NonLinear  Least  Squares  Nonlinear  Hyperbolic  Algorithm  (NLLS 
NHA).  We  have  shown  that  the  NLLS  NHA  is  a  robust  algorithm  for  constructing  OMT  fits  to 
TOVS  radiance  data,  generating  physically  meaningful  fits  on  100%  of  a  test  sample  of  45  scans. 
These  fits  were  obtained  both  on  long  wavelength  (700  cm’1)  and  short  wavelength  (2250  cm-1) 
TOVS  data. 

Acceptable  fits  to  the  TOVS  radiance  data  were  constructed  using  only  4  parameters  in 
OMT.  Such  fits  generate  OMT  temperature  profiles  with  significant  meteorological  features. 
OMT  temperature  profiles  are  most  appropriately  employed  in  the  remote  determination  of 
atmospheric  structure  parameters  and  for  the  determination  of  input  parameters  for  numerical 
weather  prediction  codes,  rather  that  replicating  the  small-scale  temperature  structures  seen  in 
radiosonde  observations. 


6.  REFERENCES 


Acton,  F.  S.  (1970):  Numerical  Methods  that  Work,  Harper  and  Row,  New  York 
Chandrasekhar,  S.  (1960):  Radiative  Transfer,  Dover,  New  York 

Hohlfeld,  R.  G.,  J.  C.  Kilian,  T.  W.  Drueding,  and  J.  F.  Ebersole  (1988):  “Novel  Methodology 
for  Application  of  Adaptive  Systems  Techniques  for  DMSP  Remote  Temperature  Sensing”, 
Creative  Optics  Report,  COl-SR-26  to  Air  Force  Geophysics  Laboratory,  Hanscom  AFB, 
AFGL-TR-88-0323 

Hopfield,  J.,  and  D.  Tank  (1985):  Biol.  Cybernetics,  52,  141 

Jeffrey,  W.  and  R.  Rosner  (1986):  “Optimization  Algorithms:  Simulated  Annealing  and  Neural 
Network  Processing”,  Astrophys.  J.,  310,  473-481 

King,  J.  I.  F.  (1985):  “Theory  and  Application  of  Differential  Inversion  to  Remote 
Temperature  Sensing”,  in  Advances  in  Remote  Temperature  Sensing  Methods,  A.  Deepak,  H. 
Fleming,  and  M.  T.  Chahine  (Eds.),  A.  Deepak  Publishing  Co.,  Williamsburg,  VA 

King,  J.  I.  F.  (1989):  “Optical  Measure  and  the  Non-Newtonian  Calculus  of  Inverse  Transfer 
Theory",  private  communication. 

King,  J.  I.  F.,  R.  G.  Hohlfeld,  and  J.  C.  Kilian  (1989):  “Application  and  Evaluation  of  a 
Differential  Inversion  Technique  for  Remote  Temperature  Sensing”,  Meteorology  and  Atmos. 
Phys.,  in  press 

Leon,  S.  J.,  and  J.  1.  F.  King  (1988):  “A  Smart  Algorithm  for  Nonlinear  Interpolation  and 
Noise  Discrimination”,  presented  at  Int’l  Workshop  on  Remote  Sensing  Retrieval  Methods, 
Williamsburg,  VA,  15-18  Dec.  ’87,  Deepak  Publ.  Co.,  in  press. 

Marquardt,  D.  W.  (1963):  J.  Soc.  Ind.  Appl.  Math.  JJ.,  pp.  431-441. 

Metropolis,  N.,  A.  Rosenbluth,  M.  Rosenbluth,  A.  Teller,  and  E.  Teller  (1953):  J.  Chem.  Phys., 
6,  1087 

Press,  W.  H.,  B.  P.  Flannery,  S.  A.  Teukolsky,  and  W.  T.  Vetterling  (1986):  Numerical  Recipes: 
The  Art  of  Scientific  Computing,  Cambridge  University  Press,  Cambridge 


27 


APPENDIX 


Remit*  of  Optical  Meuure  Theory  NLLS-NHA  applied  to  TOYS  long  wavelength  (700  cm-1)  data: 


ZONE  1 


PASS 

Parameter* 

a 

b 

K1 

LI 

Radiance 

Chi  Squared 

Temperature 

Chi  Squared 

1 

21.9897 

62.4568 

0.0712 

18.8208 

7.7253 

9.6669 

2 

20.8109 

30.8127 

0.0764 

20.7666 

8.2237 

12.6141 

3 

20.1024 

36.1872 

-0.0332 

19.7633 

10.9466 

22.1108 

4 

20.6843 

31.9178 

0.0773 

20.6341 

8.0789 

12.8196 

6 

20.2016 

42.7945 

-  0.0694 

19.6847 

8  2166 

9  0614 

e 

24.1988 

48.6561 

0.0172 

19.2390 

4.6067 

12.9135 

7 

20.0930 

41.9466 

-0.0228 

19.3296 

16.4449 

26.6887 

8 

19.7163 

42.3897 

-0.0265 

19.1709 

11.9110 

20.7660 

9 

20.2266 

46.9303 

0.0789 

19.2882 

8.6948 

7.1661 

10 

21.0669 

64.7609 

0.0634 

18.6387 

8.7642 

13.6989 

11 

21.0407 

63.2192 

0.0619 

18.7438 

8.3353 

17.3833 

12 

20.6767 

61 .3638 

-0.0291 

18.8327 

10.0276 

16.6180 

13 

20.1530 

41.7166 

-0.0187 

19.3757 

16.4502 

27.6259 

14 

20.7280 

43.4640 

0.0482 

19.6211 

6.8636 

13.6721 

16 

20.8486 

38.0441 

-0.0966 

19.7211 

7.6981 

13.1202 

ZONE  2 


PASS 

Parameter! 

a 

b 

K1 

LI 

Radiance 

Chi  Squared 

Temperature 

Chi  Squared 

i 

20.0666 

60.2691 

-0.0134 

18.4829 

21.3097 

14.0864 

2 

22.7040 

46.1479 

-0.6861 

21.0458 

26.5995 

16.1779 

3 

19.9808 

86.6364 

0.1233 

16.2895 

24.6734 

18.0511 

4 

23.6784 

91.7270 

00429 

14.7175 

47.5068 

36.3018 

6 

22.6640 

86.1693 

0.0529 

16.8039 

39.4410 

32.0148 

6 

26.4618 

92.0344 

0.0335 

13.7685 

49.9676 

40  8269 

7 

24.4512 

86.9867 

0  0361 

16.3968 

44.2754 

41  8214 

8 

23  4871 

78.6828 

0.0443 

16.7173 

35.0703 

23.7176 

9 

13.9169 

42.1279 

-0.4881 

28.6781 

13.6568 

10.4029 

10 

20.1966 

61.0701 

0.0021 

18.5280 

8.7407 

16.1631 

11 

8.8821 

33.2692 

-0.4956 

34.1463 

10.5916 

15.4403 

12 

17.8678 

40.2948 

-0.4735 

24.2218 

13.1907 

7.9640 

13 

20.1134 

76.6693 

-0.2322 

16.6688 

23.4453 

477.8820 

14 

24.7764 

91.7628 

0.0627 

14.0436 

42.4881 

33.7667 

16 

18.6641 

62.2492 

-0.4316 

23.9844 

20.4876 

12.4871 

ZONE  S 


PASS 

Parameters 

a 

b 

K1 

LI 

Radiance 

Chi  Squared 

Temperature 

Chi  Squnred 

1 

17.6488 

113.2630 

0.0060 

17.1666 

72.7358 

2 

22.6468 

113.6100 

0.0106 

13.0646 

64.9291 

3 

21.6579 

113.0420 

0.1641 

13.9088 

67.5030 

4 

20.9466 

108.1440 

-0.0603 

14.3876 

68  6499 

6 

21.1346 

110.9210 

0.0556 

14.1233 

83.5679 

6 

3.0187 

79.4961 

-0.4834 

33.9501 

80.9482 

7 

17.8404 

116.7700 

-0.0049 

16.4167 

82.3276 

8 

18.0460 

114.8840 

0.1531 

17.0912 

75.9672 

9 

18.0460 

114.8840 

0.1631 

17.0912 

75.9672 

10 

22  7688 

118.1260 

0.1193 

12.3443 

74.3194 

11 

18.6673 

112.6770 

0.0669 

16.7632 

80.9567 

12 

16.3976 

117.6640 

18  4916 

80.4202 

37.7834 

13 

22.3840 

118.3220 

12.6496 

66.9086 

66.0465 

14 

20.1937 

113.2730 

0.2993 

15.7100 

63  9672 

42.1414 

16 

22.1519 

108.6000 

-0.0683 

13.7462 

68.9993 

73.7613 

Results  of  Optical  Meaaure  Theory  NLLS-NHA  applied  to  TOVS  short  wavelength  (22SO  cm-1)  data: 


ZONE  1 


PASS 

Parameters 

a 

b 

K 1 

LI 

Radiance 

Chi  Squared 

Temperature 

Chi  Squared 

1 

-0.0817 

0.6680 

0.0136 

0.1198 

14.6169 

22.7618 

2 

-0.0318 

0.2224 

-0  0609 

0.0841 

3.62269 

16.7897 

3 

-0.0352 

0.2807 

-0.0260 

0.0856 

10.8188 

31.4143 

4 

-0.0334 

0.2322 

-0.0362 

0.0832 

4.64404 

16.6869 

6 

-0.0069 

0.3429 

0  0091 

0.0571 

9.80052 

20.2758 

6 

-0.0722 

0.6652 

0.1163 

0.1193 

16.6156 

22.4707 

7 

-0  0327 

0  3440 

-0.0110 

0.0846 

17.8616 

40.0358 

8 

-0.0391 

0.3646 

-0.0261 

0.0807 

12.8674 

32.0189 

9 

-0.0293 

0.4591 

1.1946 

0.0719 

9.6139 

14  3352 

10 

-0.0806 

0.6727 

0.1443 

0.1171 

18.3589 

26.0443 

11 

-0  0804 

0.5674 

0.1177 

0.1160 

17.2854 

26.3545 

12 

-0.0307 

0.6278 

0.0062 

0.0683 

16.6639 

28.3946 

13 

-0.0360 

0.3494 

-0.0022 

0.0823 

19.494 

38.6053 

14 

-0.0965 

0.3846 

0.0136 

0.1423 

9.51026 

24.4296 

15 

-0.2357 

0.2629 

-0.2399 

0.2944 

11.4666 

26.0404 

ZONE  2 

PASS 

Parameters 

a 

b 

K1 

_ LI 

Radiance 

Chi  Squared 

Temperature 

Chi  Squared 

1 

-0.1773 

0.4908 

-0.2828 

29.3679 

39.5405 

2 

-0.0360 

1.1118 

-46.3660 

0.0745 

7.26997 

8.8E+16 

3 

0.0230 

1.1199 

-13.1211 

0.0580 

14.9585 

61825.1 

4 

0  0983 

1.1620 

-9.8580 

0.0307 

39.7419 

3397.02 

5 

-0.0186 

1.2098 

82776200 

140509 

61.1645 

3703.7 

6 

-0.0939 

1.2368 

-2.8616 

-0.0270 

79.5673 

122.164 

7 

0.1168 

1.0662 

-9.2818 

0.0246 

32.8007 

1945.81 

8 

-0.1176 

1.1452 

3.4532 

0.2170 

82.7624 

66.6798 

9 

-0.3895 

0.4183 

-0.6167 

0.4347 

29.9299 

29.1770 

10 

-0.1438 

0.4388 

-0.4682 

0.1942 

24.3943 

43.2193 

11 

-0.6422 

0.2091 

-0.5165 

0.6007 

30.3683 

34.7671 

12 

-0.6453 

0  0244 

-0.6200 

0.7203 

32.2153 

24.5781 

13 

0.0460 

0.9901 

-11.5321 

0.0416 

12.8681 

12380.5 

14 

0.1019 

1.2926 

-3.1609 

-0.0306 

98.5834 

129.156 

16 

-0.8287 

0.2311 

-0.4774 

0.8978 

60.7268 

39.4245 

ZONE  3 
PASS 

Parameter* 

a 

b 

K1 

LI 

Radiance 

Chi  Squared 

Temperature 

Chi  Squared 

1 

-4.4988 

-0.1666 

-0.2997 

4.5656 

178.66 

120.195 

2 

0.0056 

1.8417 

-14.6069 

0.1118 

29.4934 

647087 

3 

0.0093 

1.7432 

-14.3962 

0.106B 

24.6212 

524611 

4 

0.0187 

1.8322 

-16.0639 

0.1267 

28.1216 

2196220 

6 

0.0098 

1.8338 

-16.9123 

0.1503 

27.6644 

76821500 

6 

0.0796 

1.4106 

-3.5776 

-0.0234 

161.232 

112.93 

7 

-0.0469 

1.9996 

-20.8680 

0  1666 

26.67 

2.73249E+11 

8 

-3.8966 

0.1130 

-0.2960 

3.9548 

184.67 

126.902 

9 

-3.8966 

0.1130 

-0.2960 

3.9648 

184.67 

137.040 

10 

0.0133 

1.9180 

-13.8176 

0.1147 

33.6562 

262896 

11 

-0.0260 

1.8948 

-19.3872 

0.1648 

24.927 

12826400000 

12 

-0.0762 

2.1409 

-22.3062 

0.1747 

27.8821 

4.33052E+12 

13 

0.0119 

1.9368 

-14.7067 

0  1238 

30.9684 

1080360 

14 

-4.6701 

-0.1796 

-0.2986 

4.7380 

173.639 

131.606 

16 

0.0002 

1.9937 

-16.1188 

0.1320 

32.62 

2970380 

29 


